!***********************************************************************
!                                                                      *

      SUBROUTINE DSPEVX(JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU, ABSTOL, M, W&
         , Z, LDZ, WORK, IWORK, IFAIL, INFO)
!***********************************************************************
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:04:58   1/ 3/07
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      USE cons_C
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE lsame_I
      USE dlamch_I
      USE dlansp_I
      USE dcopy_I
      USE dopgtr_I
      USE dopmtr_I
      USE dscal_I
      USE dsptrd_I
      USE dstebz_I
      USE dstein_I
      USE dsteqr_I
      USE dsterf_I
      USE dswap_I
      USE xerbla_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER  :: N
      INTEGER  :: IL
      INTEGER  :: IU
      INTEGER  :: M
      INTEGER  :: LDZ
      INTEGER  :: INFO
      REAL(DOUBLE), INTENT(IN) :: VL
      REAL(DOUBLE), INTENT(IN) :: VU
      REAL(DOUBLE), INTENT(IN) :: ABSTOL
      CHARACTER  :: JOBZ
      CHARACTER  :: RANGE
      CHARACTER  :: UPLO
      INTEGER, DIMENSION(*)  :: IWORK
      INTEGER, DIMENSION(*)  :: IFAIL
      REAL(DOUBLE), DIMENSION(*)  :: AP
      REAL(DOUBLE), DIMENSION(*)  :: W
      REAL(DOUBLE), DIMENSION(LDZ,*)  :: Z
      REAL(DOUBLE), DIMENSION(*)  :: WORK
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL, INDISP, INDIWO, &
         INDTAU, INDWRK, ISCALE, ITMP1, J, JJ, NSPLIT
      REAL(DOUBLE) :: ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, &
         SMLNUM, TMP1, VLL, VUU
      LOGICAL :: ALLEIG, INDEIG, VALEIG, WANTZ
      CHARACTER :: ORDER
!-----------------------------------------------
!   I n t r i n s i c  F u n c t i o n s
!-----------------------------------------------
      INTRINSIC MIN, SQRT
!-----------------------------------------------
!
!  -- LAPACK driver routine (version 2.0) --
!     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
!     Courant Institute, Argonne National Lab, and Rice University
!     September 30, 1994
!
!     .. Scalar Arguments ..
!     ..
!     .. Array Arguments ..
!     ..
!
!  Purpose
!  =======
!
!  DSPEVX computes selected eigenvalues and, optionally, eigenvectors
!  of a real symmetric matrix A in packed storage.  Eigenvalues/vectors
!  can be selected by specifying either a range of values or a range of
!  indices for the desired eigenvalues.
!
!  Arguments
!  =========
!
!  JOBZ    (input) CHARACTER*1
!          = 'N':  Compute eigenvalues only;
!          = 'V':  Compute eigenvalues and eigenvectors.
!
!  RANGE   (input) CHARACTER*1
!          = 'A': all eigenvalues will be found;
!          = 'V': all eigenvalues in the half-open interval (VL,VU]
!                 will be found;
!          = 'I': the IL-th through IU-th eigenvalues will be found.
!
!  UPLO    (input) CHARACTER*1
!          = 'U':  Upper triangle of A is stored;
!          = 'L':  Lower triangle of A is stored.
!
!  N       (input) INTEGER
!          The order of the matrix A.  N >= 0.
!
!  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
!          On entry, the upper or lower triangle of the symmetric matrix
!          A, packed columnwise in a linear array.  The j-th column of A
!          is stored in the array AP as follows:
!          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!
!          On exit, AP is overwritten by values generated during the
!          reduction to tridiagonal form.  If UPLO = 'U', the diagonal
!          and first superdiagonal of the tridiagonal matrix T overwrite
!          the corresponding elements of A, and if UPLO = 'L', the
!          diagonal and first subdiagonal of T overwrite the
!          corresponding elements of A.
!
!  VL      (input) DOUBLE PRECISION
!  VU      (input) DOUBLE PRECISION
!          If RANGE='V', the lower and upper bounds of the interval to
!          be searched for eigenvalues. VL < VU.
!          Not referenced if RANGE = 'A' or 'I'.
!
!  IL      (input) INTEGER
!  IU      (input) INTEGER
!          If RANGE='I', the indices (in ascending order) of the
!          smallest and largest eigenvalues to be returned.
!          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!          Not referenced if RANGE = 'A' or 'V'.
!
!  ABSTOL  (input) DOUBLE PRECISION
!          The absolute error tolerance for the eigenvalues.
!          An approximate eigenvalue is accepted as converged
!          when it is determined to lie in an interval [a,b]
!          of width less than or equal to
!
!                  ABSTOL + EPS *   max( |a|,|b| ) ,
!
!          where EPS is the machine precision.  If ABSTOL is less than
!          or equal to zero, then  EPS*|T|  will be used in its place,
!          where |T| is the 1-norm of the tridiagonal matrix obtained
!          by reducing AP to tridiagonal form.
!
!          Eigenvalues will be computed most accurately when ABSTOL is
!          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!          If this routine returns with INFO>0, indicating that some
!          eigenvectors did not converge, try setting ABSTOL to
!          2*DLAMCH('S').
!
!          See "Computing Small Singular Values of Bidiagonal Matrices
!          with Guaranteed High Relative Accuracy," by Demmel and
!          Kahan, LAPACK Working Note #3.
!
!  M       (output) INTEGER
!          The total number of eigenvalues found.  0 <= M <= N.
!          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!
!  W       (output) DOUBLE PRECISION array, dimension (N)
!          If INFO = 0, the selected eigenvalues in ascending order.
!
!  Z       (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M))
!          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!          contain the orthonormal eigenvectors of the matrix A
!          corresponding to the selected eigenvalues, with the i-th
!          column of Z holding the eigenvector associated with W(i).
!          If an eigenvector fails to converge, then that column of Z
!          contains the latest approximation to the eigenvector, and the
!          index of the eigenvector is returned in IFAIL.
!          If JOBZ = 'N', then Z is not referenced.
!          Note: the user must ensure that at least max(1,M) columns are
!          supplied in the array Z; if RANGE = 'V', the exact value of M
!          is not known in advance and an upper bound must be used.
!
!  LDZ     (input) INTEGER
!          The leading dimension of the array Z.  LDZ >= 1, and if
!          JOBZ = 'V', LDZ >= max(1,N).
!
!  WORK    (workspace) DOUBLE PRECISION array, dimension (8*N)
!
!  IWORK   (workspace) INTEGER array, dimension (5*N)
!
!  IFAIL   (output) INTEGER array, dimension (N)
!          If JOBZ = 'V', then if INFO = 0, the first M elements of
!          IFAIL are zero.  If INFO > 0, then IFAIL contains the
!          indices of the eigenvectors that failed to converge.
!          If JOBZ = 'N', then IFAIL is not referenced.
!
!  INFO    (output) INTEGER
!          = 0:  successful exit
!          < 0:  if INFO = -i, the i-th argument had an illegal value
!          > 0:  if INFO = i, then i eigenvectors failed to converge.
!                Their indices are stored in array IFAIL.
!
!  =====================================================================
!
!     .. Parameters ..
!     ..
!     .. Local Scalars ..
!     ..
!     .. External Functions ..
!     ..
!     .. External Subroutines ..
!     ..
!     .. Intrinsic Functions ..
!     ..
!     .. Executable Statements ..
!
!     Test the input parameters.
!
      WANTZ = LSAME(JOBZ,'V')
      ALLEIG = LSAME(RANGE,'A')
      VALEIG = LSAME(RANGE,'V')
      INDEIG = LSAME(RANGE,'I')
!
      INFO = 0
      IF (.NOT.(WANTZ .OR. LSAME(JOBZ,'N'))) THEN
         INFO = -1
      ELSE IF (.NOT.(ALLEIG .OR. VALEIG .OR. INDEIG)) THEN
         INFO = -2
      ELSE IF (.NOT.(LSAME(UPLO,'L') .OR. LSAME(UPLO,'U'))) THEN
         INFO = -3
      ELSE IF (N < 0) THEN
         INFO = -4
      ELSE IF (VALEIG .AND. N>0 .AND. VU<=VL) THEN
         INFO = -7
      ELSE IF (INDEIG .AND. IL<1) THEN
         INFO = -8
      ELSE IF (INDEIG .AND. (IU<MIN(N,IL) .OR. IU>N)) THEN
         INFO = -9
      ELSE IF (LDZ<1 .OR. WANTZ .AND. LDZ<N) THEN
         INFO = -14
      ENDIF
!
      IF (INFO /= 0) THEN
         CALL XERBLA ('DSPEVX', (-INFO))
         RETURN
      ENDIF
!
!     Quick return if possible
!
      M = 0
      IF (N == 0) RETURN
!
      IF (N == 1) THEN
         IF (ALLEIG .OR. INDEIG) THEN
            M = 1
            W(1) = AP(1)
         ELSE
            IF (VL<AP(1) .AND. VU>=AP(1)) THEN
               M = 1
               W(1) = AP(1)
            ENDIF
         ENDIF
         IF (WANTZ) Z(1,1) = ONE
         RETURN
      ENDIF
!
!     Get machine constants.
!
      SAFMIN = DLAMCH('Safe minimum')
      EPS = DLAMCH('Precision')
      SMLNUM = SAFMIN/EPS
      BIGNUM = ONE/SMLNUM
      RMIN = SQRT(SMLNUM)
      RMAX = MIN(SQRT(BIGNUM),ONE/SQRT(SQRT(SAFMIN)))
!
!     Scale matrix to allowable range, if necessary.
!
      ISCALE = 0
      ABSTLL = ABSTOL
      IF (VALEIG) THEN
         VLL = VL
         VUU = VU
      ENDIF
      ANRM = DLANSP('M',UPLO,N,AP,WORK)
      IF (ANRM>ZERO .AND. ANRM<RMIN) THEN
         ISCALE = 1
         SIGMA = RMIN/ANRM
      ELSE IF (ANRM > RMAX) THEN
         ISCALE = 1
         SIGMA = RMAX/ANRM
      ENDIF
      IF (ISCALE == 1) THEN
         CALL DSCAL ((N*(N + 1))/2, SIGMA, AP, 1)
         IF (ABSTOL > 0) ABSTLL = ABSTOL*SIGMA
         IF (VALEIG) THEN
            VLL = VL*SIGMA
            VUU = VU*SIGMA
         ENDIF
      ENDIF
!
!     Call DSPTRD to reduce symmetric packed matrix to tridiagonal form.
!
      INDTAU = 1
      INDE = INDTAU + N
      INDD = INDE + N
      INDWRK = INDD + N
      CALL DSPTRD (UPLO, N, AP, WORK(INDD), WORK(INDE), WORK(INDTAU), IINFO)
!
!     If all eigenvalues are desired and ABSTOL is less than or equal
!     to zero, then call DSTERF or DOPGTR and SSTEQR.  If this fails
!     for some eigenvalue, then try DSTEBZ.
!
      IF ((ALLEIG .OR. INDEIG .AND. IL==1 .AND. IU==N) .AND. ABSTOL<=ZERO) THEN
         CALL DCOPY (N, WORK(INDD), 1, W, 1)
         INDEE = INDWRK + 2*N
         IF (.NOT.WANTZ) THEN
            CALL DCOPY (N - 1, WORK(INDE), 1, WORK(INDEE), 1)
            CALL DSTERF (N, W, WORK(INDEE), INFO)
         ELSE
            CALL DOPGTR (UPLO, N, AP, WORK(INDTAU), Z, LDZ, WORK(INDWRK), IINFO&
               )
            CALL DCOPY (N - 1, WORK(INDE), 1, WORK(INDEE), 1)
            CALL DSTEQR (JOBZ, N, W, WORK(INDEE), Z, LDZ, WORK(INDWRK), INFO)
            IF (INFO == 0) THEN
               IFAIL(:N) = 0
            ENDIF
         ENDIF
         IF (INFO == 0) THEN
            M = N
            GO TO 20
         ENDIF
         INFO = 0
      ENDIF
!
!     Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
!
      IF (WANTZ) THEN
         ORDER = 'B'
      ELSE
         ORDER = 'E'
      ENDIF
      INDIBL = 1
      INDISP = INDIBL + N
      INDIWO = INDISP + N
      CALL DSTEBZ (RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL, WORK(INDD), WORK(&
         INDE), M, NSPLIT, W, IWORK(INDIBL), IWORK(INDISP), WORK(INDWRK), IWORK&
         (INDIWO), INFO)
!
      IF (WANTZ) THEN
         CALL DSTEIN (N, WORK(INDD), WORK(INDE), M, W, IWORK(INDIBL), IWORK(&
            INDISP), Z, LDZ, WORK(INDWRK), IWORK(INDIWO), IFAIL, INFO)
!
!        Apply orthogonal matrix used in reduction to tridiagonal
!        form to eigenvectors returned by DSTEIN.
!
         CALL DOPMTR ('L', UPLO, 'N', N, M, AP, WORK(INDTAU), Z, LDZ, WORK(&
            INDWRK), INFO)
      ENDIF
!
!     If matrix was scaled, then rescale eigenvalues appropriately.
!
   20 CONTINUE
      IF (ISCALE == 1) THEN
         IF (INFO == 0) THEN
            IMAX = M
         ELSE
            IMAX = INFO - 1
         ENDIF
         CALL DSCAL (IMAX, ONE/SIGMA, W, 1)
      ENDIF
!
!     If eigenvalues are not in order, then sort them, along with
!     eigenvectors.
!
      IF (WANTZ) THEN
         DO J = 1, M - 1
            I = 0
            TMP1 = W(J)
            DO JJ = J + 1, M
               IF (W(JJ) >= TMP1) CYCLE
               I = JJ
               TMP1 = W(JJ)
            END DO
!
            IF (I == 0) CYCLE
            ITMP1 = IWORK(INDIBL+I-1)
            W(I) = W(J)
            IWORK(INDIBL+I-1) = IWORK(INDIBL+J-1)
            W(J) = TMP1
            IWORK(INDIBL+J-1) = ITMP1
            CALL DSWAP (N, Z(1,I), 1, Z(1,J), 1)
            IF (INFO == 0) CYCLE
            ITMP1 = IFAIL(I)
            IFAIL(I) = IFAIL(J)
            IFAIL(J) = ITMP1
         END DO
      ENDIF
!
      RETURN
!
!     End of DSPEVX
!
      END SUBROUTINE DSPEVX
